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ABSTRACT 


A  new  unstructured  grid  two-dimensional,  depth- integrated  (2DDI),  morphodynamic 
model  is  presented  for  the  prediction  of  morphological  evolutions  in  shallow  water.  This 
modelling  system  consists  of  two  coupled  model  components:  i.)  a  well  verified  and  validated 
continuous  Galerkin  (CG)  finite  element  hydrodynamic  model;  and  ii.)  a  new  sediment 
transport/bed  evolution  model  that  uses  a  discontinuous  Galerkin  (DG)  method  for  the  solution  of 
the  sediment  continuity  equation.  The  DG  method  is  a  robust  finite  element  method  that  is 
particularly  well  suited  for  this  type  of  advection  dominated  transport  equation.  It  incoiporates 
upwinded  numerical  fluxes  and  slope  limiters  to  provide  sharp  resolution  of  steep  bathymetric 
gradients  that  may  form  in  the  solution,  and  it  possesses  a  local  conservation  property  that 
conserves  sediment  mass  on  an  elemental  level.  In  this  paper,  we  focus  specifically  on  the 
implementation  and  verification  of  the  DG  model.  Details  are  given  on  the  implementation  of  the 
method,  and  numerical  results  are  presented  for  three  idealized  test  cases  which  demonstrate  the 
accuracy  and  robustness  of  the  method  and  its  applicability  in  predicting  medium-term 
morphological  changes  in  channels  and  coastal  inlets. 


1. 


INTRODUCTION 


The  transport  of  sediment  as  bed  load  is  an  important  process  that  occurs  in  rivers, 
estuaries,  and  coastal  regions.  In  many  situations,  this  process  and  the  resulting  morphological 
changes  of  the  bed  can  have  a  detrimental  impact  on  the  coastal  infrastructure  and  environment. 
For  example,  dredged  navigational  channels  and  coastal  inlets  can  be  rendered  almost  entirely 
useless  by  the  accumulation  of  transported  sediment.  Returning  these  structures  to  operational 
status,  through  dredging  operations  or  the  construction  of  jetties,  represents  a  significant  cost  to 
the  agencies  that  maintain  them.  As  another  example,  the  structural  integrity  of  bridges  and  piers 
may  be  compromised  due  to  excessive  scour  of  the  bed  around  abutments.  In  addition  to  these 
infrastructure  problems,  there  is  host  of  environmental  issues  of  concern,  such  as  the  transport  of 
pollutants  with  or  as  sediment,  that  can  cause  serious  ecological  damage.  Accurate  numerical 
models  that  can  predict  sediment  transport  and  the  resulting  bed  morphology  can  help  manage 
these  costly  problems. 

Clearly,  the  processes  of  sediment  transport  and  morphological  evolution  of  the  bed  are 
determined  by  the  properties  of  the  fluid  flow,  which  in  turn  are  affected  by  the  changes  in  the 
morphology  of  the  bed  that  they  induce.  Thus,  the  motion  of  the  fluid  and  the  motion  of  the  bed 
form  an  interdependent  two-phase  phenomenon  that  must  be  analyzed  using  a  model  system  made 
up  of  two  distinct  but  interdependent  model  components:  i.)  a  hydrodynamic  component  defining 
the  evolution  of  the  flow,  and  ii.)  a  sediment  transport/morphological  component  defining  the 
evolution  of  the  bed.  Such  a  modelling  system  is  often  referred  to  as  a  moiphodynamic  model.  A 
description  and  comparison  of  some  existing  morphodynamic  model  systems  is  given  by 
Nicholson,  et  al.  (1997).  Typically,  these  model  systems  use  structured  computational  grid 
methods.  To  a  lesser  extent,  unstructured  grid  methods  have  also  been  implemented  and  can,  in 
fact,  be  highly  advantageous  based  on  their  ability  to  provide  local  grid  refinement  near  important 
bathymetric  features  and  structures.  The  ability  to  provide  local  grid  refinement  where  it  is 
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needed  leads  to  improved  accuracy  for  a  given  computational  cost  as  compared  to  models  that  use 
structured  grid  methods.  However,  both  structured  and  unstructured  grid  method  solutions  to  the 
governing  morphological  equation  can  experience  numerical  robustness  and  accuracy  problems 
manifested  in  the  form  of  spurious  spatial  oscillations,  especially  in  the  presence  of  steep 
bathymetric  gradients  (see  for  example  Johnson  and  Zyserman,  2002). 

In  this  paper,  we  describe  the  development  of  a  new  unstructured  grid  morphodynamic 
model  system  that  uses  a  new  class  of  highly  accurate  finite  element  methods  for  the  solution  of 
the  governing  morphological  equation.  The  hydrodynamic  model  component  of  our  system  is 
provided  by  the  well  verified  and  validated  unstructured  grid  model  ADCIRC,  developed  by  the 
second  author  and  a  number  of  collaborators  (Luettich  and  Westerink,  2004).  ADCIRC  is  both  a 
two-dimensional,  depth-integrated  (2DDI)  and  three-dimensional  (3D)  free  surface  flow  model. 
In  this  paper,  we  focus  specifically  on  the  2DDI  ADCIRC  model,  which  solves  the  shallow  water 
equations  using  the  standard  or  continuous  Galerkin  (CG)  finite  element  method  in  space.  To 
overcome  well  known  problems  in  solving  the  shallow  water  equations  using  equal-order 
interpolating  spaces  with  the  CG  finite  element  method,  the  continuity  equation  is  replaced  by  the 
so-called  generalized  wave  continuity  equation  (GWCE)  (Lynch  and  Gray,  1979  and  Kinnmark, 
1986).  The  solution  strategy  used  in  ADCIRC  has  proven  to  be  robust  and  computationally 
efficient,  and  it  has  been  validated  in  a  large  number  of  cases  (see  for  example  Blain,  et  al.,  1994; 
Westerink  et  al.,  1994;  Mukai,  et  al.,  2002;  Westerink  et  al.,  2004). 

Working  with  a  well  established  hydrodynamic  model  then,  the  main  focus  of  this  paper 
is  the  development  and  verification  of  a  bed  load  sediment  transport/morphological  model 
component  to  work  in  conjunction  with  ADCIRC.  Mathematically,  the  morphological  evolution 
of  the  bed  is  defined  by  the  so-called  sediment  continuity  or  Exner  equation.  This  equation 
simply  states  that  the  time  rate  of  change  of  the  bed  elevation  is  equal  to  the  divergence  of  the 
sediment  flux,  which  can  be  expressed  in  terms  of  the  local  flow  properties  through  the  use  of  an 
empirical  sediment  transport  formula.  As  is  well  known,  solving  advection  dominated  transport 
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equations  of  this  type  using  the  CG  finite  element  method  will  frequently  lead  to  spurious  spatial 
oscillations  in  the  solution.  To  overcome  these  shortcomings,  a  number  of  so-called  advection 
schemes  can  be  employed  (see  Iskandarani,  et  al.,  2005  for  a  review  and  comparison  of  some  of 
the  more  popular  schemes).  One  such  scheme  that  has  received  considerable  recent  attention  and 
that  has  been  applied  successfully  to  a  wide  variety  of  problems  is  the  discontinuous  Galerkin 
(DG)  finite  element  method. 

Originally  developed  by  Reed  and  Hill  (1973),  but  more  recently  expounded  on  in  a 
series  of  papers  by  Cockbum,  et  al.  (see  the  review  article  by  Cockbum  and  Shu,  2001  and  the 
references  therein),  the  DG  method  uses  trial  and  test  function  spaces  that  are  continuous  over  a 
given  element  but  which  allow  discontinuities  between  elements.  This  results  in  a  block  diagonal 
or,  with  an  appropriate  choice  of  basis,  diagonal  mass  matrix  that  can  be  trivially  inverted. 
Communication  between  elements  is  accomplished  via  a  so-called  numerical  flux,  which  for  the 
case  of  a  scalar  equation  can  be  defined  using  upwinding  techniques.  The  method  is  also  “locally 
conservative”,  meaning  that  the  conservation  of  the  transported  quantity  is  satisfied  on  a  local  or 
elemental  level.  This  has  been  shown  to  be  a  desirable  property  when  coupling  flow  and 
transport  algorithms  (see  for  example  Dawson,  et  al.,  2004) 

In  this  paper,  we  present  the  implementation  and  verification  of  a  DG  sediment 
transport/morphological  model  that  is  coupled  to  the  ADCIRC  hydrodynamic  model.  We  note 
that  this  sediment  transport  model  is  just  one  component  of  a  suite  of  DG  model  components  that 
are  currently  being  developed  for  flow  and  transport,  which  will  form  a  completely  DG  based 
moiphodynamic  modeling  system  with  both  h  (grid  size)  and  p  (polynomial  order)  refinement 
options.  In  this  paper,  we  restrict  our  attention  to  the  second-order  (p  =  1)  case  for  the  sediment 
transport  model,  but  we  note  that  /> refinement  is  easily  implemented  within  the  framework  of  the 
DG  method  (see  Kubatko,  et  al.,  2005  for  an  example  of  this  for  the  shallow  water  equations). 

This  paper  is  organized  as  follows.  In  Section  2,  we  describe  the  mathematical  model 
defining  the  sediment  transport  and  morphological  evolution  of  the  bed  which  consists  of  the 
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sediment  continuity  equation  and  an  empirical  sediment  transport  formula.  We  then  present  a 
simplified  mathematical  model,  which  we  refer  to  as  the  Exner  model,  that  uncouples  the 
sediment  transport/morphological  model  from  the  hydrodynamic  model.  This  simplified  model 
can  be  used  as  a  verification  tool  for  the  numerical  method.  In  Section  3,  we  give  a  detailed 
description  of  our  implementation  of  the  DG  method  for  the  sediment  continuity  equation,  giving 
specific  details  on  the  numerical  flux,  basis,  quadrature  rules,  time  discretization,  slope  limiter, 
and  continuous  projection  that  are  employed.  In  Section  4,  we  present  numerical  results  from 
three  test  cases  with  the  aim  of:  i.)  verifying  that  the  method  achieves  second-order  convergence 
in  space,  and  ii.)  demonstrating  how  the  model  can  be  used  for  predicting  so  called  medium-term 
(see  for  example  de  Vriend,  et  al.,  1993)  morphological  changes  in  channels  and  coastal  inlets. 
Finally,  in  Section  5,  we  summarize  this  paper,  and  we  briefly  discuss  the  current  and  future  work 
in  the  development  of  this  model  system. 
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2. 


MATHEMATICAL  MODEL 


The  evolution  of  the  bed  or  bottom  surface  elevation  due  to  the  transport  of  sediment  as 
bed  load  is  governed  by  the  so-called  sediment  continuity  or  Exner  equation  (see  for  example 
Henderson,  1966): 


-  +  v.qi  =  0  (1) 
at 


where  z  is  the  elevation  of  the  bed  relative  to  a  datum  located  below  the  bed  (z  is  positive 
upwards)  and  qA  is  the  bed  load  sediment  transport  function  vector. 

In  order  to  close  Eq.  (1),  a  functional  form  of  q*  must  be  specified.  It  is  assumed  that  the 
sediment  transport  is  always  in  the  direction  of  the  flow  velocity,  U  =  ( u ,  v)  where  u  and  v  are  the 
velocity  components  of  the  flow  in  the  x  and  y  directions,  respectively.  Thus  the  vector  qA  is 
computed  as: 


q*  =U|q6|  (2) 

where  |q4 1  is  the  magnitude  of  the  sediment  transport  in  the  direction  of  the  flow  and  U  is  the 

normalized  flow  velocity  vector  (i.e.  U  =  U/|U|).  There  are  a  number  of  empirical  bed  load 

sediment  transport  functions  available  (e.g.  Bagnold,  Einstein,  Meyer-Peter  and  Mueller,  see  for 
example  Sleath,  1984  for  a  thorough  list),  most  of  which  can  be  transcribed  in  the  following  form: 

|q;,|  =  H(U,//,---)|U|"  (3) 
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where  A  is  a  given  function  and  n  is  a  given  positive  constant  both  of  which  are  specific  to  the 
particular  sediment  transport  formula.  Note  that  A  is  typically  a  function  of  the  flow  velocity,  U, 
the  total  height  of  the  water  column,  H  =  'Q  -  z  (where  C  is  the  water  surface  elevation  relative  to 
the  same  datum  as  the  bed),  and  a  number  of  constants  that  are  based  on  sediment  properties  such 
as  sediment  type  and  grain  size  and  data  fitting  procedures.  The  constant  n  is  typically  in  the 
range  of  1  <  n  <  3  . 

In  our  model,  we  will  make  use  of  a  new  bed  load  formula  developed  by  Camenen  and 
Larson  (2005),  though  the  numerical  model  to  be  described  will  be  general  enough  to  allow  the 
use  of  any  sediment  transport  formula  provided  it  is  a  function  of  H  and  a  monotonically 
increasing  function  of  U.  Camenen  and  Larson  develop  new  bed  load  sediment  transport 
formulas  for  transport  due  to  currents,  waves,  and  combined  waves  and  currents.  Their  formulas 
were  shown  to  provide  the  best  agreement  with  the  data  sets  that  were  compiled  compared  to  a 
number  of  previously  proposed  formulas  (Camenen  and  Larson,  2005). 

In  this  paper,  we  only  consider  the  Camenen  and  Larson  bed  load  sediment  transport 
formula  due  to  currents  which  is  given  by  (in  dimensional  form  -  SI  units): 


1  1  „  i  5 

f  T  \ 

m  =  Czc  -  exp 

-4.5  — 

l  C  J 

(4) 


where  zc  is  the  shear  stress  at  the  bottom  due  to  the  current,  zcr  is  the  critical  shear  stress,  and  C  is 
a  constant  given  by: 


C  = 


12 

gy[p{p,  -  P ) 


(5) 
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where  g  is  the  acceleration  due  to  gravity  and  p  and  ps  are  the  water  and  sediment  density, 
respectively.  The  shear  stress  is  computed  by  the  formula: 


r 


c 


(6) 


where  /  is  a  dimensionless  friction  factor  calculated  assuming  a  logarithmic  velocity  profile  (see 
for  example  Sleath,  1984): 


/(") 


(  d<n  V 

1  +  In 

50 

1.15/7 

(7) 


where  dsn  is  the  median  grain  size.  The  critical  shear  stress  is  computed  from  the  critical  Shields 
parameter  which  is  either  estimated  as  a  constant  or  calculated  using  the  formula  proposed  by 
Soulsby  and  Whitehouse  (1997). 

Using  Eqs.  (5)  -  (7)  the  sediment  transport  formula  can  be  written  in  the  form  of  Equation 
(3)  with  n  =  3  and  A  given  by  the  function: 


A  =  C 


(1  5 

1.5 

f  \ 

r 

-Pf 

exp 

-4.5  — 

U  J 

l  Tc) 

(8) 


Note  that  A  is  a  monotonically  increasing  function  of  rc  and  therefore  U. 
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Rigid  Lid 


u 

l=>  H=z-Z 


1 

2  I -  Bed  Datum 

_ J _  _ 

Fig.  1:  Definition  sketch  for  Exner’s  model 


2.1  A  simplified  model 

For  purposes  of  verifying  our  numerical  scheme,  we  use  a  simplified  mathematical  model 
that  essentially  uncouples  the  sediment  continuity  equation  from  the  hydrodynamics.  This  allows 
us  to  verify  the  underlying  numerics  of  the  model  in  a  simplified  setting  by  comparing  it  to 
analytical  solutions.  Assume  that  the  flow  is  unidirectional  (say  in  the  x  directional  only)  and 
quasi-steady  with  a  rigid  lid.  With  these  assumptions,  the  flow  velocity  is  given  by: 


fi_  qf 

~  H  ~  fi-z 


(9) 


where  qj  is  a  constant  flow  discharge  and  £  is  the  elevation  of  the  rigid  lid  measured  from  the  bed 
datum  (see  Fig.  1).  Furthermore,  assume  that  the  sediment  transport  is  given  by  Eq.  (3)  with  A  = 
constant  and  n  =  1.  With  these  assumptions  the  sediment  continuity  equation  can  be  written  as: 


dz  d 

- 1 - 

dt  dx 


r 


v 


=  0 


(10) 
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This  model  was  originally  proposed  by  Exner  (1925).  Assuming  a  smooth  initial  condition, 


z(x,  0)  =  z0 ,  the  classical  solution  is  given  implicitly  by: 


z(x,t)  =  z0(x-cj),  c.  =  Aqf/[<Z -zf  (11) 

where  cz  is  the  propagation  speed  of  the  bed.  As  is  well  known,  non-linear  hyperbolic  equations 
such  as  Eq.  (10),  depending  on  the  initial  conditions,  will  develop  steep  gradients  (and  eventually 
discontinuities  or  shocks)  which  provide  a  rigorous  test  for  a  numerical  method.  A  similar  model 
was  examined  by  Johnson  and  Zyserman  (2002)  in  the  context  of  testing  a  finite  difference 
scheme. 
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Fig.  2:  A  typical  element  e  and  its  neighboring  element  along  edge  i  with  normal  ni;  v""'  and  vexj  denote 
the  value  of  a  function  v  along  edge  i  when  approaching  the  edge  from  the  interior  and  exterior  of  the 
element  respectively. 


3.  NUMERICAL  MODEL 

In  this  section,  we  give  a  detailed  description  of  our  DG  sediment 
transport/morphological  model.  To  begin  we  define  some  notation.  Given  a  spatial  domain,  G, 
which  has  been  discretized  into  a  set  of  non-overlapping  elements,  let  D.e  define  the  domain  of  a 
typical  element  e  and  denote  the  boundary  of  the  element  by  T(, .  Our  numerical  approximation  of 
z  will  make  use  of  piecewise  smooth  Junctions  which  are  continuous  over  G(  but  which  allow 
discontinuities  between  elements  along  a  given  edge.  We  denote  this  space  of  functions  by  Vh. 

Given  a  smooth  function  v  defined  over  e,  we  denote  the  values  of  v  along  an  edge  by  v*"'1  when 

approaching  the  edge  from  the  interior  of  the  element  and  v^when  approaching  the  edge  from 
the  exterior  of  the  element.  The  outward  unit  normal  vector  for  the  boundary  of  the  element  will 
be  denoted  by  n,  and  the  fixed  unit  normal  vector  for  a  given  edge  i  will  be  denoted  by  n,  (see 
Fig.  2). 

In  our  numerical  scheme,  we  will  also  make  use  of  continuous,  piecewise  linear 
approximations  of  U  and  C  obtained  from  the  ADCIRC  model  to  compute  the  local  sediment 
transport  rates.  Briefly,  these  approximations  are  obtained  by  solving  the  shallow  water 
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equations  using  the  CG  Galerkin  finite  element  method  in  space  and  implicit/explicit  time 
stepping  (see  Luettich  and  Westerink  (2004)  for  details).  As  previously  mentioned,  to  achieve 
non-oscillatory  results  the  primitive  continuity  equation  is  replaced  with  the  GWCE. 

We  apply  the  DG  method  to  the  sediment  continuity  equation  by  multiplying  Eq.  (1)  by 
a  test  function  v  e  Vh  and  integrating  over  Qe  to  obtain: 


(12) 


Integrating  the  second  term  of  this  equation  by  parts  gives: 


L^vdne  -  L Vv  ■ q,//Q  +  L v  q>  ■ m/r  =  0 


(13) 


Next  we  replace  the  solution  z  with  an  approximate  solution  zh  which,  using  Galerkin’ s 
method,  is  constructed  from  a  set  of  basis  functions  which  belong  to  the  same  space,  Vh,  as  the 
test  functions.  Due  to  the  fact  that  there  may  be  discontinuities  along  element  edges,  the 
boundary  integral  of  Eq.  (13)  is  undefined  and  for  this  we  define  a  numerical  flux,  qh .  In  our 
formulation,  we  use  a  simple  upwind  flux  based  on  the  assumption  that  the  sediment  transport  is 
in  the  direction  of  the  current: 


jq6W-n,  U-n  >  0 
I  q*<ev) ' n,  U-n  <  0 


(14) 


With  the  approximate  solution  and  the  numerical  flux  defined,  the  weak  formulation  of  the 
problem  now  becomes: 
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b  lH  Vd^e  _  VV '  +  Jr. V  & dF «  =  ° 


(15) 


Note  that  the  method  is  locally  or  elementally  conservative  in  the  following  sense:  setting  v  =  1 
on  £!,,  and  zero  elsewhere  we  have: 


Li7dn-+ U‘dr-=0  (16) 

That  is,  the  time  rate  of  change  of  zh  over  Q,  is  balanced  by  the  net  flux  of  sediment  into  Qe. 

We  proceed  by  describing  the  details  of  the  implementation  of  the  scheme  including  the 
choice  of  basis  functions,  the  quadrature  rules  employed  to  compute  the  integrals,  the  time 
discretization,  the  application  of  a  slope  limiter  to  eliminate  local  undershoots  or  overshoots  in 
the  solution  in  the  presence  of  steep  gradients,  and  the  continuous  projection  procedure  used  to 
project  the  discontinuous  approximation  zh  into  the  space  of  continuous,  piecewise  linear 
functions  which  are  fed  back  into  ADCIRC  as  updated  bathymetry. 
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h, 


h2 

Fig.  3:  Master  element  defined  in  local  coordinates  <f  and  t]  showing  the  degrees  of  freedom  /?,  -  the 
value  ofh  at  the  midpoint  of  edge  i  opposite  of  corner  node  i. 

3. 1  Basis  and  degrees  of  freedom 

As  emphasized  by  Cockbum  and  Shu  (1998),  we  note  here  that  a  judicious  choice  of 
basis  functions  can  simplify  the  implementation  of  the  scheme  and  improve  the  computational 
efficiency.  Owing  to  the  fact  that  discontinuities  are  permitted  across  element  interfaces,  the 
choice  of  the  basis  functions  are  not  limited  by  the  requirement  of  continuity  as  in  the  CG  finite 
element  method.  Therefore,  one  can  choose  degrees  of  freedom  that,  for  example,  save  cost  in 
evaluating  the  integrals  in  Eq.  (15)  and/or  simplify  the  implementation  of  the  slope  limiter.  In  our 
implementation,  we  use  piecewise  linear  triangular  elements  described  below. 

Considering  the  “master  element”  as  shown  in  Fig.  3  defined  in  the  transformed 
coordinates  f  and  17,  the  approximate  solution  zh  can  be  expressed  as: 

07) 

1=1 
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where  the  degrees  of  freedom,  z,  are  the  values  of  the  approximate  solution  at  the  mid-point  of 
each  edge  and  the  basis  functions,  fat  define  the  linear  element  of  Crouzeix  and  Raviart  (1973) 
which  for  the  master  element  shown  in  Fig.  3  can  be  written  in  the  form: 

fa=\-24,  fa2-\-  2ij,  fa  -2^  +  2p  -  \  (18) 

There  are  several  things  to  note  about  this  basis.  The  functions  fa  are  equal  to  1  at  the  mid-point 
of  each  edge  i  and  0  at  the  mid-points  of  the  other  two  edges.  The  basis  functions  are  orthogonal 
over  an  element,  specifically: 


i*j 

where  Qm  denotes  the  domain  of  the  master  element.  This  property,  of  course,  gives  rise  to  an 
orthogonal  mass  matrix  that  can  be  trivially  inverted.  Lastly,  in  the  continuous  projection 
procedure  to  be  described,  we  will  make  use  of  the  value  of  the  approximate  solution  at  the 
vertices  of  the  triangle.  The  value  of  zh  at  vertex  i,  denoted  by  zvi,  which  is  the  vertex  opposite  of 
edge  i  (see  Fig.  3),  is  easily  computed  as: 


Jn  fafajd^dii  = 


[1/6, 

lo, 


z«=-2zi+YJzJ  (20) 

j= 1 

As  a  final  note,  we  remark  that  the  orthogonal,  hierarchical,  “modal”  type  basis  proposed 
by  Dubiner  (1991),  which  simplifies  p  refinement  and  also  adaptivity,  can  easily  be  implemented 
within  the  framework  of  the  DG  method. 
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3.2  Quadrature  Rules 


Both  of  the  integrals  appearing  in  Eq.  (15)  are  evaluated  using  suitable  numerical 
quadrature  rules.  We  note  that  by  using  numerical  quadrature  and  the  simple  upwind  numerical 
flux  defined  previously,  we  can  easily  implement  a  number  of  different  sediment  transport 
formulas  into  the  scheme  without  making  any  changes  to  the  base  algorithm  itself  (provided  that 
the  formula  meets  the  requirements  as  specified  in  Section  2).  Cockbum  and  Shu  (1998)  note 
that  for  a  DG  spatial  discretization  of  degree  p ,  quadrature  rules  that  are  exact  for  polynomials  of 
degree  2 p  and  2p+\  must  be  used  for  the  area  and  boundary  integrals,  respectively.  Thus  for  the 
linear  elements  used  here  (p  =  1)  we  use  a  three  point  quadrature  rule  for  the  triangle  so  the  area 
integral  of  Equation  (15)  is  approximated  by  (noting  that  Vv  is  constant  over  the  element): 


Vv 


\D  Q  A )  *  Vv •  (  £  w,.q/;  (U „ ) 


V  1=1 


(21) 


where  the  w,’s  are  the  quadrature  weights  of  the  associated  quadrature  points,  which  are  the 
midpoints  of  each  edge.  Using  this  rule,  the  sediment  transport  function,  qA  is  easily  evaluated  at 
the  quadrature  points  given  the  fact  that  we  already  have  z„  which  are  the  degrees  of  freedom,  and 
we  need  only  to  compute  U  and  C  at  the  mid-point  of  each  edge.  We  note  that  these  values  are 
easily  obtained  by  averaging  the  two  vertices  for  the  given  edge  (owing  to  the  fact  that  U  and  C 
are  approximated  using  linear  functions  over  the  element  as  well).  The  boundary  integrals,  which 
must  integrate  a  third  degree  polynomial  exactly,  are  evaluated  using  the  two-point  Legendre- 
Gauss  quadrature  rule. 
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3.3 


Time  Discretization 


The  DG  spatial  discretization  reduces  the  problem  to  a  system  of  ordinary  differential 
equations  which  we  write  in  the  concise  form: 


jt{z„)  =  Lk{zh,\Sh&h)  (22) 

where  zh  is  the  vector  of  unknowns  over  the  whole  domain. 

We  discretize  this  system  of  equations  in  time  using  a  second-order  Runge-Kutta  scheme, 
which  is  equivalent  to  the  so-called  modified  Euler  method,  written  in  the  form: 


z(1)=z(,)  + 

*-h  L h  ^ 


,(,+1)  _ 


AtmLh(zy,\jU,&) 

=  i(zW  +  zf»+A«.J,(zf'.oi'>,Cl")) 


(23) 


where  A tm  is  the  moiphological  time  step  which  may  be  different  than  that  of  the  hydrodynamic 
time  step,  A th ,  and  where  it  is  to  be  noted  that  U  and  C,  are  held  fixed  at  time  t. 

Given  that  explicit  time  stepping  is  used,  the  size  of  the  morphological  time  step  is 
limited  by  a  Courant-Friedrichs-Levy  (CFL)  condition.  A  direct  calculation  of  this  condition 
proves  difficult  in  practice  due  to  the  highly  non-linear  nature  of  the  sediment  transport  function, 
and  instead  we  simply  take  A tm  -  N  x  Ath ,  where  N  is  some  positive  integer  usually  in  the  range 
of  10  to  50,  i.e.  the  bed  is  updated  every  10  to  50  hydrodynamic  time  steps.  In  practice,  this 
approach  has  proven  to  work  well  for  a  wide  variety  of  problems  and  requires  little  additional 
computational  effort.  It  has  been  estimated  that  using  this  approach  the  additional  computational 
cost  for  running  the  morphodynamic  model  is  on  the  order  of  2-10%  of  the  cost  of  running  the 
hydrodynamic  model  alone. 
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3.4  Slope  Limiter- 


In  order  to  prevent  spurious  oscillations  at  sharp  fronts,  a  slope  limiter  is  applied  at  each 
step  of  the  Runge-Kutta  method  described  above.  We  apply  a  simple  slope  limiter  in  which  the 
degrees  of  freedom  z,  for  a  given  element  e  are  compared  to  the  average  of  the  approximate 
solution  over  e,  zeavg  and  the  average  of  the  neighboring  element  e'  of  the  given  edge,  zeavg  .  If  z, 

does  not  fall  in  between  the  values  zeavg  and  zeavg  for  the  given  edge  i  then  the  degrees  of  freedom 

for  element  e  are  set  equal  to  zeavg  .  In  this  way,  the  average  of  the  element  is  maintained  while 

setting  its  slope  equal  to  zero,  and  sediment  mass  is  still  conserved  over  the  element.  It  should  be 
noted  this  slope  limiter  is  very  easy  to  implement,  but  it  can  cause  some  numerical  smoothing  of 
the  solution.  More  sophisticated  limiters  that  are  less  dissipative  are  currently  being  investigated. 

We  remark  that  for  sufficiently  smooth  bathymetries,  in  practice  it  is  often  unnecessary  to 
apply  the  limiter.  However,  as  the  bed  evolves,  steep  gradients  may  develop  in  the  bed,  and  it  has 
been  observed  that  without  the  use  of  a  limiter  oscillations  develop  in  the  neighborhood  of  the 
steep  gradient.  Typically,  however,  these  oscillations  seem  to  remain  localized  and  do  not 
degrade  the  solution  globally.  The  role  of  the  slope  limiter  then,  at  least  for  the  problems 
examined,  is  that  of  a  mechanism  to  eliminate  local  oscillations  rather  than  for  stabilizing  the 
scheme. 
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Fig.  4:  Nodey  surrounded  by  n  elements;  the  continuous  approximation  z.  tor  node  /  is  determined 
from  the  (possibly)  n  unique  values  from  the  elements  surrounding  the  node. 


3.5  Continuous  projection  procedure 


As  previously  mentioned,  ADCIRC  makes  use  of  approximations  that  are  continuous  in 
space  across  the  entire  domain.  Thus,  in  addition  to  our  discontinuous  approximation  Z/„  we  must 
also  compute  a  continuous  approximation  which  must  be  fed  back  into  ADCIRC  after  computing 
the  updated  bathymetry.  We  wish  to  accomplish  the  following  with  our  procedure:  given  a  node 
j  which  is  a  vertex  for  n  different  elements  we  wish  to  compute  a  single  nodal  value  denoted  by 
z .  based  on  the  n  (possibly)  unique  values  at  that  node  that  are  obtained  from  the  DG  method 

within  the  individual  elements  attached  to  that  node  (see  Fig.  4).  We  have  experimented  with 
several  different  approaches  for  obtaining  these  single  nodal  values  and  based  on  numerical 
experiments  have  implemented  an  angle  based  weighted  average  given  by: 


0=Z 


— A,- 

V  '-SUM  ) 


(24) 


where  Z;.  is  the  angle  of  the  vertex  of  element  i  and  ZS(/M  is  the  total  sum  of  the  vertex  angles 
around  node  j.  We  have  also  experimented  with  weighted  area  averaging  and  centroidal  type 
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averaging,  but  we  have  found  that  the  approach  given  by  Equation  (24)  gives  the  most  consistent 
results  under  a  wide  variety  of  grid  configurations.  We  note  that  under  certain  grid  configurations 
it  was  observed  that  mild  in-plane  (x-y  plane)  oscillations  appeared  in  the  continuous 
representation  of  the  bed.  The  weighted  angle  approach  minimized  the  appearance  of  these 
oscillations,  which  were  often  much  more  visible  using  other  averaging  techniques.  It  should  be 
noted  this  procedure  does  not  affect  the  local  conservation  property  of  the  sediment  due  to  the 
fact  that  z  is  not  actually  used  in  the  computations  for  updating  the  bed. 
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4. 


NUMERICAL  RESULTS 


The  DG  method  outlined  above  has  been  applied  to  a  number  of  problems.  In  this 
section,  we  show  the  results  for  three  idealized  test  cases. 

4.1  Test  Casel:  Morphological  evolution  of  a  symmetric  mound 

In  this  test  case,  we  apply  the  DG  method  outlined  above  to  the  Exner  model  introduced 
in  section  2.  The  Exner  model  is  examined  in  order  to  verify  the  numerical  method 
independently  of  the  hydrodynamic  model.  It  also  affords  us  the  opportunity  to  compare  our 
numerical  results  to  exact  solutions  so  that  we  may  check  the  order  of  convergence  of  the  method. 

We  solve  a  problem  originally  posed  by  Exner  (1925).  The  problem  examines  the 
evolution  of  an  initially  symmetric  mound  subjected  to  steady,  uni-directional  flow  with  a  rigid- 
lid  assumption  for  the  flow.  The  initial  condition  is  given  as: 


:  (x,  y,  0)  =  z0  (x,  y)  =  A0  +  A1  cos 


2  nx 


(25) 


where  the  parameters  A0  ,  Aj  ,  and  2  are  as  defined  in  Fig.  5  which  shows  a  cross  section  of  the 
mound  along  the  x-axis.  We  take  A0  =  A1  =  1,  2  =  20  in  Equation  (25)  and  if  =  3,  Aqj  =  1  in  Eq. 
(10).  The  flow  is  assumed  to  be  in  the  x  direction  only,  and  we  use  periodic  boundary  conditions, 
i.e.  z(x  =  -/t/2,  v)  =  z(x  =  +X/2,  v)  and  z[x,y  =  -zl/2)  =  z(x,y  =  +A/2)  The  exact 
solution  is  given  by  Eq.  (11). 
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*  =  -A/2 


_  A 
x  =  0 


x  =  A/2 


Fig.  5:  Cross  section  of  the  initial  condition  for  Exner’s  “dune  problem”  as  defined  by  Eq.  (25) 


We  solve  this  one-dimesional  problem  over  a  two-dimensional  domain  using  four 
different  grids  with  uniform  nodal  spacing  of  h  =  1,  0.5,  0.25,  and  0.125.  We  compute  the 
maximum  or  L&  norm  by  comparing  our  DG  solutions  to  the  exact  solution.  In  Fig.  6,  we  plot  h 
versus  the  maximum  error  norm  on  a  log-log  scale  where  it  can  be  observed  that  the  theoretical 
convergence  rate  of  p+ 1  is  obtained.  Both  the  numerical  and  exact  solutions  of  the  evolution  of 
the  mound  at  a  cross  section  taken  along  the  x-axis  are  shown  at  two  different  times  in  Fig.  7. 
The  solutions  indicate  that  the  mound  develops  into  a  dune-like  shape  with  a  gentle  upstream 
slope  followed  by  a  steeper  downstream  slope  that  becomes  progressively  steeper  in  time.  It 
should  be  noted  how  well  the  DG  solution  captures  the  steep  downstream  slope  of  the  dune 
without  the  introduction  of  any  spurious  spatial  oscillations  or  any  significant  numerical  damping. 
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Bed  Elevation  Convergence 


Fig.  6:  Convergence  plot  of  test  case  1  demonstrating  2nd  order  convergence 


Fig.  7:  Comparison  of  the  exact  and  DG  solution  for  test  case  1 
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B  =  250  m 


Fig.  8:  Computational  grid  of  test  case  2 


4.2  Test  Case  2:  Converging  channel 

In  this  problem,  we  return  to  the  full  motphodynamic  modelling  system  and  examine  the 
morphological  evolution  of  an  initially  flat  bed  in  a  converging  channel.  A  plan  view  of  the 
channel  showing  the  computation  grid  is  shown  in  Fig.  8.  The  channel  tapers  in  from  a  maximum 
width  of  500  m  at  the  edges  to  250  m  in  the  center  over  a  distance  of  2  km.  The  boundary 
conditions  for  the  hydrodynamics  are  specified  in  such  a  way  that  a  maximum  velocity  of 
approximately  1  m/s  occurs  in  the  center  of  the  channel.  The  evolution  of  the  bed  is  examined 
over  a  90  day  period.  The  sediment  density  and  median  grain  size  of  the  bed  are  taken  to  be  2000 
kg/nT  and  0.2  mm,  respectively.  The  time  step  used  in  the  hydrodynamic  model  is  2  seconds  and 
the  bed  is  updated  every  50  hydrodynamic  time  steps.  Figs.  9a-9c  show  plots  of  the  bed  elevation 
surface  and  velocity  contours  at  30,  60,  and  90  days.  The  bed  changes  have  been  scaled  in  the 
vertical  for  easy  visualization. 

The  velocity  throughout  the  channel  varies  from  approximately  0.50  m/s  at  the  ends  of 
the  channel  to  approximately  1  m/s  in  the  center  of  the  channel.  The  bed  experiences  erosion  in 
the  converging  part  of  the  channel  due  to  the  increase  in  the  flow  velocity.  Conversely,  in  the 


23 


diverging  part  of  the  channel,  as  the  flow  velocity  decreases,  accretion  of  the  sediment  occurs  and 
a  mound  or  shoal  develops  in  time.  It  can  be  noted  the  scour  and  accretion  patterns  occurring  in 
the  center  of  the  channel  are  slightly  larger  than  those  occurring  toward  the  sides  of  the  channel 
across  the  width  of  a  given  cross  section.  This  can  be  explained  by  the  fact  that  the  velocity  field 
is  not  entirely  uniform  across  the  width  of  the  channel  with  somewhat  higher  velocities  occurring 
in  the  center.  These  small  variations  in  the  velocity  field  across  the  width  of  the  channel  produce 
variations  in  the  morphology  of  the  bed  across  the  width  of  the  channel  given  the  fact  that  the 
sediment  transport  is  a  function  of  U3.  We  also  note  that  the  velocity  field  evolves  along  with  the 
morphological  changes. 

Finally,  we  remark  that  the  computed  results  of  the  evolution  of  the  bed  compare  well 
qualitatively  to  an  analytical  solution  given  by  Exner  (1925)  for  a  problem  of  the  same  geometry. 
Exner’s  results,  as  shown  in  Figure  10,  are  the  solution  of  a  simplified  model  similar  to  that  of 
Section  2  but  modified  accordingly  to  account  for  variations  in  the  width  of  the  channel  (see  Graf, 
1971  for  details).  Specifically,  it  can  be  seen  that  the  numerical  and  analytical  solutions  show  the 
same  general  evolution  of  the  bed,  i.e.  scour  in  the  converging  section  of  the  channel  and 
accretion  in  the  diverging  section. 
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Fig.  9a:  Velocity  contours  and  bed  surface  for  test  case  2  -  Day  30 
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Fig.  9b:  Velocity  contours  and  bed  surface  for  test  case  2  -  Day  60 
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Fig.  9c:  Velocity  contours  and  bed  surface  for  test  case  2  -  Day  90 
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28 


Fig.  11:  a.)  Computational  grid  for  the  idealized  inlet  of  test  case  3  and  b.)  details  of  the  grid  in  the 
vicinity  of  the  inlet. 

4.3  Test  Case3:  Idealized  Inlet 

In  this  problem,  we  apply  the  model  to  the  case  of  an  idealized  inlet  as  shown  in  Fig. 
1 1 .  The  domain  consists  of  a  1 0  km  by  20  km  sound  connected  to  the  open  ocean  through  an 
inlet  which  is  1  km  wide  and  0.5  km  long.  The  open  ocean  boundary  is  20  km  from  the  entrance 
of  the  inlet  and  is  50  km  wide.  The  initial  bathymetry  in  the  sound  and  through  the  inlet  is 
constant  at  a  depth  of  5  m.  South  of  the  inlet  the  bathymetry  varies  linearly  from  5  m  at  the 
entrance  to  14  m  at  the  open  ocean  boundary.  The  sediment  density  and  median  grain  size  are  the 
same  as  those  specified  in  the  previous  problem.  The  grid  for  the  problem  is  also  shown  in  Fig. 
1 1  with  the  inset  showing  the  details  in  the  vicinity  of  the  inlet.  The  nodal  spacing  ranges  from 
100  m  near  the  inlet  to  1  km  at  the  open  ocean  boundary.  The  problem  is  forced  with  an  M2  tide 
with  a  15  cm  amplitude  which  produces  a  maximum  current  of  approximately  1  m/s  through  the 
inlet.  The  time  step  used  for  the  hydrodynamics  is  5  seconds,  and  the  bed  is  updated  every  50 
hydrodynamic  time  steps.  A  28  day  simulation  was  run  with  magnified  sediment  transport  rates 
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in  order  to  enhance  the  advective  processes  and  accelerate  the  morphological  evolution  of  the 
dominant  features.  Runs  with  no  magnification  of  the  sediment  transport  rates  produced 
qualitatively  similar  results  with  smaller  changes  in  the  bed. 

Figs.  12a-12d  show  the  time  evolution  of  the  bed  every  7  days  over  the  28  day 
simulation.  Note  that  larger  values  in  the  bed  elevation  indicate  erosion  and  lower  values  indicate 
accretion  due  to  the  fact  that  the  bed  is  measured  as  positive  downward  from  the  geoid.  On  day 
7,  there  is  noticeable  erosion  beginning  at  the  southern  end  of  the  inlet.  Accumulation  of  the 
sediment  can  be  seen  along  the  sides  of  the  inlet  and  to  the  south  of  the  inlet  indicating  the  initial 
formation  of  an  ebb  shoal.  During  flood  tide  on  day  14,  it  can  be  seen  that  there  has  been 
significant  erosion  through  the  throat  of  the  inlet  resulting  in  the  initial  formation  of  a  flood  shoal. 
It  can  also  be  seen  that  the  ebb  shoal  has  become  more  pronounced.  By  day  21,  there  are  distinct 
flood  and  ebb  shoals  to  the  north  and  south  of  the  inlet  respectively.  There  is  also  additional 
erosion  through  the  inlet  following  the  same  pattern  as  the  initial  scour.  At  the  end  of  day  28, 
there  has  been  significant  scour  through  the  entire  length  of  the  inlet  and  the  flood  and  ebb  shoals 
have  become  even  more  pronounced.  It  should  be  noted  that  even  at  this  level  of  coarse  grid 
resolution  the  model  captures  the  main  morphological  changes  one  expects  to  observe  in  tidally 
dominated  coastal  inlets  (see  for  example  Hayes,  1980). 
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Fig.  12.)  Evolution  of  the  bed  in  the  vicinity  of  the  inlet  over  the  28  day  simulation  -  a.)  Day  7, 

b.)  Day  14,  c.)  Day  21,  and  d.)  Day  28 
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5. 


SUMMARY  AND  FUTURE  WORK 


In  this  paper,  we  have  presented  a  new  unstructured  grid  morphodynamic  model  which 
makes  use  of  the  existing  ADCIRC  finite  element  hydrodynamic  model  and  a  new  DG  finite 
element  sediment  transport/morphological  model.  Specific  details  were  given  on  the 
implementation  of  the  DG  method,  and  the  model  was  shown  to  produce  good  results  in  three 
idealized  test  cases.  In  the  first  test  case  it  was  verified,  through  the  use  of  the  Exner  model,  that 
the  method  achieves  second-order  convergence  in  space.  Additionally,  it  was  demonstrated  how 
the  DG  method  can  accurately  capture  steep  gradients  in  the  bathymetry  without  the  introduction 
of  spurious  spatial  oscillations.  The  second  and  third  test  cases  demonstrated  how  the  full 
moiphodynamic  modelling  system  can  be  used  to  predict  medium-term  morphological  changes  of 
the  bed  in  channels  and  tidally  dominated  coastal  inlets. 

We  conclude  with  some  comments  on  the  current  development  of  this  morphodynamic 
modelling  system,  in  terms  of  both  physical  and  numerical  features  that  will  be  implemented.  In 
this  paper,  we  have  only  considered  sediment  transport  due  to  currents.  However,  in  many 
coastal  scenarios  short  waves,  which  interact  with  the  current  through  the  introduction  of 
radiation  stress  terms  in  the  momentum  equations,  can  be  the  dominant  force  in  the  sediment 
transport  process.  Therefore,  future  work  will  involve  coupling  a  wave  model  component  into  the 
modelling  system  to  include  the  effects  of  waves  in  both  the  hydrodynamics  and  sediment 
transport  processes.  Numerically,  as  was  previously  indicated,  the  present  model  is  only  one 
component  of  a  suite  of  DG  models  that  are  currently  being  developed.  Other  DG  model 
components  will  include  a  2DDI  DG  hydrodynamic  model  (see  Kubatko,  et  al.,  2005)  and  2DDI 
DG  transport  models  for  salinity  and  temperature.  In  many  applications,  these  models  will  be 
used  in  advection  dominated  flow  scenarios  such  as  coastal  inlets.  The  DG  method  is  particularly 
advantageous  for  these  types  of  situations. 
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